function  [out,errorMes]  = biasCorrectTheta2(theta2Hat,Var_u,Cov10_u,Cov01_u,factors,...
    bandwidthT,allowDeSaling,numBootStep2,seedNum,numBootStep1And3)
% This function bias-corrects the estimates of Theta2 using a simple bootstrap procedure.
dispOn           = 1;
nonStatScalingOn = 1;

% Resampling
rng(seedNum,'twister')
[nx,T]   = size(factors);
% resampling from 2,3,...,T. That is we exclude the first observation
% because here residuals, cov10, cov01 are missing
%indexBoot         = randi(T-1,numBootStep2,T)+1;

% Unfold theta2Hat
[h0Hat,hxHat] = unfoldTheta2(theta2Hat,nx);
residuals     = zeros(nx,T);
for t=2:T
    residuals(:,t) =  factors(:,t) - (h0Hat + hxHat*factors(:,t-1));
end
theta2HatValues     = struc2values(theta2Hat);
theta2BootValues    = zeros(size(theta2HatValues,1),numBootStep2);
AVarTheta2BootDraws = zeros(size(theta2HatValues,1),size(theta2HatValues,1),numBootStep1And3);

count = 0;
countNonStat = 0;
for s=1:numBootStep2
    % resampling from 2,3,...,T. That is we exclude the first observation
    % because here residuals, cov10, cov01 are missing
    indexBoot         = randi(T-1,1,T)+1;

    % Generating the data for bootstap s
    factors_Boot     = zeros(nx,T);
    Var_u_Boot       = zeros(nx,nx,T);
    Cov10_u_Boot     = zeros(nx,nx,T);
    Cov01_u_Boot     = zeros(nx,nx,T);
    factors_Boot(:,1)= factors(:,indexBoot(1,1));
    Var_u_Boot(:,:,1)= Var_u(:,:,indexBoot(1,1));
    for t=2:T
        factors_Boot(:,t)   = h0Hat + hxHat*factors_Boot(:,t-1) + residuals(:,indexBoot(1,t));
        Var_u_Boot(:,:,t)   = Var_u(:,:,indexBoot(1,t));
        Cov10_u_Boot(:,:,t) = Cov10_u(:,:,indexBoot(1,t));
        Cov01_u_Boot(:,:,t) = Cov01_u(:,:,indexBoot(1,t));
    end
    % Using the estimator on the bootstrap sample
    [theta2Boot,~,check,~,hxBoot]  = Var1_GMM_closedForm(Var_u_Boot,Cov10_u_Boot,Cov01_u_Boot,factors_Boot,allowDeSaling);
    if check == 0
        count = count + 1;
        theta2BootValues(:,count) = struc2values(theta2Boot);
        if any(eig(hxBoot) >= 1)
            countNonStat = countNonStat + 1;
        end
    end
    
    if s <= numBootStep1And3 
        % We compute the co-variance matrix of theta2Boot for the bootstrap in step 3 of the SR approach
        AVarTheta2Step2Boot = getAVarTheta2(theta2Boot,Var_u_Boot,Cov10_u_Boot,Cov01_u_Boot,factors_Boot,bandwidthT);
        AVarTheta2BootDraws(:,:,s) = AVarTheta2Step2Boot.VarRobust;
    end
end

% The bias adjustment - ensuring stationarity of the adjusted estimate
theta2BootAvgValues = mean(theta2BootValues,2);
theta2BootAvg       = values2struct(theta2BootAvgValues,fieldnames(theta2Hat));
biasTheta2Values    = theta2BootAvgValues-theta2HatValues;
biasTheta2          = values2struct(biasTheta2Values,fieldnames(theta2Hat));
%This corresponds to: theta2BiasAdjValues = 2*theta2HatValues-theta2BootAvgValues;
theta2AdjValues     = theta2HatValues-biasTheta2Values; 
theta2Adj           = values2struct(theta2AdjValues,fieldnames(theta2Hat));

[~,hxAdj]           = unfoldTheta2(theta2Adj,nx);
eigValues           = eig(hxAdj);
realEigValues       = real(eigValues);
imagEigValues       = imag(eigValues);
modulus             = sqrt(realEigValues.^2 + imagEigValues.^2);
errorMes            = 0;
if nonStatScalingOn == 1
    if any(modulus >= 1)
        [h0AdjScaled,hxAdjScaled,sigmaAdjScaled,scalingOpt] = runInduceStatVAR(hxAdj,factors,Var_u,Cov10_u,Cov01_u);
        theta2Adj = setTheta2(h0AdjScaled,hxAdjScaled,sigmaAdjScaled);
        if dispOn == 1
            disp(['The BIAS-adjustement is scaled by ', num2str(scalingOpt),' to induce stationarity']);
            disp(['Unscaled eigenvalues = ',num2str(modulus')]);
            disp(['Scaled eigenvalues   = ',num2str(eig(hxAdjScaled)')]);
        end
    else
        % No scaling is needed, but we need to update h0 and sigma
        scalingOpt     = 1;
        meanFactors    = mean(factors,2);
        input.varXData = (factors-repmat(meanFactors,1,T))*(factors-repmat(meanFactors,1,T))'/(T-1)-mean(Var_u,3);
        input.Var_u    = Var_u;
        input.Cov10_u  = Cov10_u;
        input.Cov01_u  = Cov01_u;
        input.factors  = factors;
        input.hx       = hxAdj;
        [~,h0Adj,hxAdj,sigmaAdj] = induceStatVAR(scalingOpt,input);
        theta2Adj = setTheta2(h0Adj,hxAdj,sigmaAdj);
    end   
else
    if any(modulus >= 1)
        errorMes = 1;
        if dispOn == 1
            disp('The VAR model is nonstationary with the bootstrap adjustment');
            disp(['Eigenvalues = ',num2str(eigValues')]);
        end
    end
end

% Co-variance matrix of theta2 from bootstrap
VarBoot = (repmat(theta2BootAvgValues,1,numBootStep2)-theta2BootValues)...
    *(repmat(theta2BootAvgValues,1,numBootStep2)-theta2BootValues)'/numBootStep2;

% Re-centering the draws to account for the bias-adjustment
theta2BootDraws = (theta2BootValues - repmat(theta2BootAvgValues,1,numBootStep2)) ...
    + repmat(struct2array(theta2Adj)',1,numBootStep2);

% The t-tests (under homoskedasticity)
numTheta2          = size(theta2BootValues,1);
tValuesTheta2Asym  = zeros(numBootStep1And3,numTheta2);
for k=1:numBootStep1And3
    tValuesTheta2Asym(k,:)  = (theta2BootDraws(:,k)' - struct2array(theta2Adj))./sqrt(diag(AVarTheta2BootDraws(:,:,k)))';
end
prctilesValues          = [1.0 2.5 5.0 10 50 90 95 97.5 99];
theta2_prctiles         = zeros(numTheta2,size(prctilesValues,2)); 
tTestTheta2_prctiles    = zeros(numTheta2,size(prctilesValues,2)); 
AbstTestTheta2_prctiles = zeros(numTheta2,size(prctilesValues,2)); 
for i=1:numTheta2
    theta2_prctiles(i,:)         = prctile(theta2BootDraws(i,:)',prctilesValues);
    tTestTheta2_prctiles(i,:)    = prctile(tValuesTheta2Asym(:,i),prctilesValues);
    AbstTestTheta2_prctiles(i,:) = prctile(abs(tValuesTheta2Asym(:,i)),prctilesValues);
end
% Storing the percentiles in a struct
for j=1:size(prctilesValues,2)
    label = num2str(prctilesValues(1,j));
    label = strrep(label,'.','_');
    theta2_pTile.(['pct',label]) = theta2_prctiles(:,j);
    tTestTheta2_pTile.(['pct',label])    = tTestTheta2_prctiles(:,j);
    AbstTestTheta2_pTile.(['pct',label]) = AbstTestTheta2_prctiles(:,j);
end
   
%Saving output
out.theta2Hat               = theta2Hat;
out.theta2BootAvg           = theta2BootAvg;
out.biasTheta2              = biasTheta2;
out.theta2Adj               = theta2Adj;
out.eigValues               = eigValues;
out.countNonStat            = countNonStat;
out.VarBoot                 = VarBoot;
out.theta2BootDraws         = theta2BootDraws;
out.AVarTheta2BootDraws     = AVarTheta2BootDraws;
out.prctilesValues          = prctilesValues;
out.theta2_pTile            = theta2_pTile;
out.tTestTheta2_pTile       = tTestTheta2_pTile;
out.AbstTestTheta2_pTile    = AbstTestTheta2_pTile;
end

